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Abstract. Two counter-propagating cool and equally dense electron beams are 
modelled with particle-in-cell (PIC) simulations. The electron beam filamentation 
instability is examined in one spatial dimension. The box length resolves one pair of 
current filaments. A small, a medium-sized and a large filament are considered and 
compared. The magnetic field amplitude at the saturation time of the filamentation 
instability is proportional to the filament size. It is demonstrated, that the force on the 
electrons imposed by the electrostatic field, which develops during the nonlinear stage 
of the instability, oscillates around a mean value that equals the magnetic pressure 
gradient force. The forces acting on the electrons due to the electrostatic and the 
magnetic field have a similar strength. The electrostatic field reduces the confining 
force close to the stable equilibrium of each filament and increases it farther away. The 
confining potential is not sinusoidal, as assumed by the magnetic trapping model, and 
it permits an overlap of current filaments (plasmons) with an opposite flow direction. 
The scaling of the saturation amplitude of the magnetic field with the filament size 
observed here thus differs from that expected from the magnetic trapping model. The 
latter nevertheless gives a good estimate for the magnetic saturation amplitude. The 
increase of the peak electrostatic and magnetic field amplitudes with the filament size 
implies, that the electrons heat up more and that the spatial modulation of their mean 
speed along the beam flow direction increases with the filament size. 
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1. Introduction 

The electron beam filamentation instability (FI) generates magnetic fields in energetic 
astrophysical and solar flare plasmas [H [21 El HI E] • The FI is also relevant for inertial 
confinement fusion (ICF) [HI 0, [8] and laboratory astrophysics [9] experiments, as well 
as for particle accelerators pUl E]. The FI is driven by counter-propagating electron 
beams and it feeds on their mean flow energy. This contrasts the Weibel instability, 
which grows magnetic fields at the expense of a thermal anisotropy [T2l [T3l [Hj. The 
Weibel instability and the FI can be combined to form cumulative instabilities [15] . The 
FI becomes important, if the beam speeds are at least mildly relativistic. 

The growth and saturation of the FI can be modelled with particle-in-cell (PIC) or 
Vlasov codes. The FI has been investigated with a one-dimensional (ID) Vlasov code 
[T71 [T8] and with a two-dimensional (2D) PIC code [19], taking into account the ion 
dynamics. The impact of a flow-aligned magnetic field has been examined in Ref. [20J. 
PIC simulation studies have addressed the statistical properties of the FI in ID [21 J and 
in 2D [22l [23] . The counterstreaming electron instability has also been examined with 
3D PIC simulations [24]. 

The probably simplest and thus widely researched plasma configuration that gives 
rise to the FI consists of equally dense and equally warm electron beams that have a 
Maxwellian velocity distribution. Their thermal velocity spread Vth in the rest frame of 
the respective beam is the same in all directions. This system can be considered in a 
simulation reference frame, in which both beams move into opposite directions at the 
non-relativistic speed modulus i>& and with v t h <C v b . The FI can be isolated by selecting 
a ID or 2D simulation box, that is oriented orthogonally to the beam velocity vector v&. 
The electron velocities must be resolved in the simulation direction or plane and along 
the beam direction. The FI competes in reality with the two-stream modes and it can 
probably not be observed experimentally in an isolated form, even if the equal beam 
densities favor the FI over the two-stream instability [25]. However, the gained insight 
into the development of the isolated FI will improve the understanding of systems, in 
which the FI interplays with other instabilities. 

The linear and nonlinear evolution of the FI driven by counter-propagating identical 
electron beams is qualitatively understood, at least in one spatial dimension where the 
filament mergers are suppressed [16]. The Refs. [HI [161 [HI I2Q1 [21] have provided an 
insight into its linear and nonlinear development, which can be summarized as follows. 
The FI triggers the aperiodic growth of waves out of an initial perturbation (noise) with 
the wavevectors k _L v& over a wide band oik = |k|, where the wavenumbers k are of the 
order of the inverse electron skin depth. The electrons are deflected by the magnetic field 
perturbation, and electrons moving in opposite directions separate in space. The net 
current of these flow channels amplifies the initial perturbation and, thus, the tendency 
to form current channels. The magnetic field amplitude grows exponentially. Magnetic 
trapping has been identified as a possible saturation mechanism |14| . It has also been 
proposed [18J that the electrostatic fields, which grow during the nonlinear evolution of 
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the FI, may be important for the saturation. These electrostatic fields have been related 
to the magnetic pressure gradient [20j EE] • However, it has not yet been demonstrated 
quantitatively that the electrostatic fields during the quasi-linear evolution of the FI 
do originate from the magnetic pressure gradient. A direct comparison of the relative 
importance of the electric and magnetic fields for the nonlinear saturation of the FI is 
also lacking and this is an objective of this paper. 

The length of the ID simulation box with periodic boundaries can be selected such, 
that only one wave period grows. This FI evolves like that in a much larger box [21] 
and the saturation mechanism should thus be the same. We employ here a simulation 
box that resolves a single spatial period of the growing wave and we can thus compare 
our results to previous work [18]. We employ plasma parameters that are identical to 
those in Ref. [21] and resort to the distribution of the filament sizes, which is computed 
in that paper. We perform three simulations, in which we vary the box length. The 
spectrum of unstable wavenumbers of the FI is bounded at low and large wavenumbers, 
the latter by thermal effects |18j . The bounded /c-spectrum implies in turn a maximum 
and a minimum filament size. We model a filament size close to the minimum value, 
one close to the maximum value and one, that is close to the size with the maximum 
probability. We compare the properties of the filaments. 

This paper is structured as follows. Section 2 discusses the PIC code and the initial 
conditions. The results are presented in the section 3. All simulations demonstrate that 
the electrons are heated up orthogonally to v^, in line with previous simulations [21J. 
The heating is achieved by the simultaneous interaction of the electrons with the quasi- 
static magnetic field and the oscillatory electrostatic fields. The heating is much stronger 
for the larger filaments than for the small one. The magnetic field amplitude reaches a 
value that is in reasonable agreement with the one expected from the magnetic trapping 
mechanism. It does, however, not scale correctly with the filament size. A reason is 
that the electrostatic field modifies the potential. The force excerted by the electrostatic 
field in the simulation of the least turbulent small filament oscillates around a mean 
value that equals the magnetic pressure gradient force, confirming experimentally their 
connection. The fields of the two larger filaments show the same spatial correlation. 
The heated electrons cannot be confined by the electromagnetic field but a cooler, 
dense electron population remains localized, forming a plasmon. This plasmon can 
propagate at a sizeable fraction of the initial electron thermal speed, which contrasts 
the non-propagating filamentation modes out of which the plasmon forms. The plasmon 
maintains the net current along and, thus, the magnetic field. The slow extraordinary 
mode is pumped in the two large simulation boxes. However, the resulting growth of 
the discrete wave spectrum observed here, which has been reported first by Ref. |18j . is 
a finite box effect. The results are discussed in more detail in section 4. 
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2. Solved equations and initial conditions 

The particle-in-cell simulation method is detailed in Ref. [26]. Our code is based on 
the numerical scheme proposed by [27]. A PIC code approximates a phase space fluid 
by an ensemble of computational particles (CPs). The CPs can have a mass m cp and 
charge q cp that differs from the physical particle it represents, but the charge-to-mass 
ratio must be preserved. The equations, which the PIC code is solving, are 

dB dE 

V x E = -— , V x B = fj, J + (1) 

V-E = p/e , V-B = 0, (2) 

-jjVcp = Qcp (E[xpp] + v cp x B[x c J) , — x cp = v cp , (3) 

with p cp = m cp r c pV C p. The currents j C p = q cp v cp of each CP are interpolated to the grid. 
The summation over all CPs gives the macroscopic current J, which is defined on the 
grid. The J updates the E, B fields through Eq.(l). The updated fields are interpolated 
to the position of each CP and update its position x cp and p cp through Eq.(3). This 
scheme is repeated over N s time increments At. 

Two equally dense beams of electrons with q cp /m cp = —e/m e move along the z- 
direction. Beam 1 has the mean speed v&i = i>&z and the beam 2 has = — v &z with 
vt/c = 0.3. Both beams have a Maxwellian velocity distribution in their respective rest 
frame with a thermal speed v t h = {hT /m e ) ' 5 of Vf,/v t h = 18. Both beams are spatially 
uniform. The ID simulation box with its periodic boundary conditions is aligned with 
the x-direction. We thus denote positions by the scalar x. The plasma frequency of each 
beam with the number density n e is u p = (e 2 n e /m e eo)°' 5 . The total plasma frequency 
fip = \2u) p . The electric and magnetic fields are normalized to E^ = eE/ 'cm e Vt p and 
Bjv = eT5/m e Q p and the current to J at = 3/2n e ec. The physical position and time are 
normalized as xn = x/\ s with the electron skin depth X s = c/Q p and — tQ p . We 
drop the indices N and x, t, E, B, J are specified in normalized units. 

The x-direction is resolved in the three simulations by N g = 500 grid cells with 
the length A x . The phase space distributions fi(x, p) of beam 1 and f2(x,p) of beam 
2 are each sampled by N p = 6.05 • 10 7 CPs. The total phase space density is defined 
as f(x,p) = f\(x,p) + f 2 {x,p). The box length of run 1 is L\ — 2.8, that of run 2 
is L 2 = 2 and that of run 3 is L 3 = 0.89. We use Li,L 2 ,L 3 to label the simulation 
runs. According to the size distribution in Ref. [21], we may expect that a single wave 
period grows in each box. The simulation L 3 captures the smallest and L\ the largest 
filament that can grow with a significant probability. The simulation L 2 corresponds 
to a filament size close to that, which grows with the maximum probability. Multiple 
small filaments of a size ~ L 3 could grow in the larger simulation boxes. This is initially 
observed in the largest box L%, but the smaller filaments merge to give one large one. 
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b) E / Time c) E / Time 

Figure 1. Field amplitudes: (a) shows the B y (ki,t) for the box L3 (thick solid line), 
the box L2 (thin solid line) and the box L\ (dashed line). The amplitudes grow 
exponentially, they saturate and oscillate around an equilibrium value. This value 
increases with the box size, (b) shows E x (k2,t) for the box L2 and (c) that for the 
box L3. The E x component evolves towards an equilibrium value for the box L\ (not 
shown) and L2 and E x remains oscillatory for the box L3. 



3. Simulation results 

The selected beam velocity vector || z and the simulation box orientation imply, that 
the electrons of both beams and their micro-currents are re-distributed by the FI only 
along x. The charge- and current-neutral plasma is transformed into one with J z (x) 7^ 0. 
According to Eq.(l) a growth of J z (x) is coupled only into the E z (x) and B y (x) field 
components. This is, because the gradients along the y, z-direction vanish in the ID 
geometry. Ampere's law simplifies to dB y /dx = J z + dE z /dt. A J z oc sin {kx) gives a 
E z oc sin {kx) and B y oc — cos (/ex) so that E z and B y will have a phase shift of 90°. 
The electron re-distribution leads to a space charge and thus to a E x (x) 7^ 0. 

The B y (x,t) and E x (x,t) are Fourier transformed over space by B y (k n ,t) = 
Ng l Y^j=\ B y (jA x , t) exp (— ijk n A x ) with k n = 2nn/N g A x . The E x (x,t) is transformed 
accordingly. Figure [1] displays the time-evolution of the amplitude moduli of the 
dominant fci-mode for B y (k,t) and of the /c2-Hiode for E x (k,t). The B y (ki,t) grows 
exponentially at the rate J\i = fij ^ — 0.25 Q p until t ~ 40 in Li, L 2 and at ^3 = 0.22 Q p 
until t w 50 in the L3. The maximum linear growth rate obtained from a cold fluid model 
is fl it f = fl p Vb/cT{vb) ~ 0.29 Q p [20]. The measured growth rate is reduced compared to 
Qij in particular for the simulation L 3 by thermal effects, which cause damping at large 
k. The amplitude of B y (ki,t) saturates and remains almost constant in all simulations 
after t = 80. The saturation amplitudes are B y (k\,t > 80, Li) rs 0.06, B y (k\,t > 
80, L2) ~ 0.045 and B y (ki,t > 80, L3) w 0.02. The magnetic trapping mechanism sets 

1 /2 

in, when Q is j is comparable to the magnetic bounce frequency Qb — (ekiVbB y / 'm e ) 1 in 
physical units [14]. The Qb for the measured growth rates above are Qb/Q p = 0.2 for 
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Figure 2. The field amplitudes in the box L\\ The panels (a-c) show B y , E x and E z , 
respectively. The B y - amplitude reaches a steady state value, which convects slowly in 
space. The E x and E z components are oscillatory in space and time. The E x oscillates 
in space twice as fast as B y and the maxima of both are co-moving. The E z shows 
a phase shift of 90° compared to B y when the fields saturate. The B y ,E x fields at 
t — 50 are displayed by (d) and they show a correlation only at x w 1.3. 



Li and 0.21 for L 2 and L 3 . The agreement is excellent for the simulation L 3 , but the 
increase of Qij in the simulations Li,L 2 compared to that in L 3 does not change Qb- 
What we find instead is, that B y (ki, t > 80, Zy) oc Lj. 

The E x (k 2 ,t) grows when B y (ki,t) has reached a large amplitude and the growth 
rate of E x (k 2 ,t) is twice that of B y (ki,t), as it has previously been reported [HJ I2~T] . 
The E x (k 2 , t) in the simulation L 3 oscillates between a peak value and zero, whereas we 
find damped oscillations around a steady state value in L 2 ,L\ after t ~ 50. The peak 
amplitude of E x in the simulation L 2 exceeds that of L 3 by the factor 3 and that in Li 
is even stronger (not shown). 

3.1. Simulation L\ 

The evolution of the electromagnetic fields and of the related electron phase space is now 
examined in more detail. Figure [2] shows the relevant fields in space and time for the 
simulation L\. It also compares the spatial distributions of B y (x,t) and E x (x,t) at the 
time t = 50, when the FI saturates. The simulation box fits one oscillation of B y (x,t) 
at any fixed time t after the saturation, but the more rapid oscillations along x during 
40 < t < 50 indicate that at least initially several modes are competing. Eventually, 
the current channels merge and form a steady state distribution in ID [HJ [21]. The 
magnetic field structure slowly convects to larger x. The phase speed of this structure 
is constant after t = 70 and it amounts to ~ Vth/5. The oscillations of E x (x, t) are more 
rapid than those of B y (x,t) for any fixed 50 < t < 55. The maxima of E x (x,t) after 
t = 70 show a clear correlation with the structure of B y (x,t), because both convect 
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Figure 3. The 10-log of the phase space densities in units of CPs at t = 50 (a-c) 
and t = 120 (d-f) in the box L\. Panels (a,d) show the f(x,p z ) with p = m e VbT(vb). 
The beam temperature along p z is unchanged and the < p z {x) >\.2 are obviously 
modulated in space. The density oscillates by the factor w 10 2 . The fi(x,p x ) is shown 
in (b,e) and f2(x,p x ) in (c,f). The electrons of both beams spatially separate and (e,f) 
reveal dense turbulent electron clouds immersed in a tenuous hot electron background, 
which reaches a thermal spread « po ■ 

with the same phase speed. The amplitudes of B y (x,t = 50) and E x (x,t = 50) show 
a possible correlation only at x ~ 1.3. The smaller filaments are merging to the larger 
one at this time, complicating the interpretation of the field correlation. The E z (x,t) 
has the same spatial oscillation frequency as B y (x, t) for any fixed 40 < t < 55 and their 
phase is shifted by the expected 90°. The amplitude of E z is significantly lower than 
that of E x . The oscillations of E z (x,t) after t = 70 show no pronounced structures. 

Figure [3] displays the phase space densities fi t 2(x,p x ) and f(x,p z ) at the times 
t = 50 and t = 120. They show a significant modulation and the density maxima of 
both beams are shifted by L±/2 at t — 120. The values of fi^i^iPx) an d f(x,p z ) vary 
by two orders of magnitude. The mean values < p z {x) >i,2= / Pzfi,2(x, p) dp for each 
beam show the same oscillatory modulation with x as in Ref. [21J. The f{x,p z ) at t = 50 
also shows that several filaments develop. For example, a density maximum is found at 
p x ~ 0.2 and p z = po that is spatially correlated with a minimum at p z = —po- The 
Po = m e VbT(vb). The absolute minima of f2(x,p x ) at x ~ 1.5 and of fi(x,p x ) at x ~ 2.5 
are also not shifted by L\/2. The smaller filaments cause the rapid spatial modulation 
of B y (x,t = 50) in Fig. [21 The electrons are heated up in p x by the saturation of the 
FI from an initial thermal spread of p x /po ~ 0.05 to a peak value of p x ~ po- The 
summation of f(p x ) = J f(x, p) dxdp y dp z over many filaments will give a distribution, 
that decreases exponentially over a wide range of p x [21j. 

The supplementary movie 1 shows the 10-logarithmic f\(x,p x ,t) and f\(x,p z ,t) in 
the simulation L\. It demonstrates that only the core electrons in Fig. [3]remain spatially 
confined. This dense core population maintains B y (x,t > 50) in Fig. [2j The structures 



Saturation of the filamentation instability 



8 




12 1 2 

c) E -field / Position d) Position 



Figure 4. The field amplitudes in the box L^: The panels (a-c) show B y , E x and E z , 
respectively. The By-amplitude reaches a steady state value, which convects slowly 
in space. The E x and E z components are oscillatory in space and time. The E x 
oscillates in space twice as fast as B y , which is obvious at t w 50. The maxima of 
both are co-moving. The E z shows a phase shift of 90° compared to B y when the 
fields saturate. The B y , E x fields at t — 50 are displayed by (d) and both show a clear 
spatial correlation. The E x is modulated by its first harmonic. 



in f(x,p x ) in the movie 1, one of which occurs in Fig. [3fe) at x ~ 0.7 and p x ~ 0, 
resemble phase space holes [28]. The potentials of these structures contribute to the 
E x (x,t), causing its rapid fluctuations in Fig. [2H). These phase space holes complicate 
the identification of the relation between E x , B y in Fig. EH) in addition to the ongoing 
merging of small filaments to larger ones. 

3.2. Simulation L 2 

We reduce now the box length from L\ = 2.8A S to L 2 = 2X S . The comparison of 
Fig. H] with Fig. [2] reveals some similarities between the respective field evolutions. 
The structures in the saturated B y (x,t) and E x (x,t) fields are co- moving also in the 
simulation L 2 and they have here a phase speed ~ —v t h/G- The E z (x,t) and B y (x,t) 
have the same spatial oscillation period when they saturate at t w 40, again shifted 
by 90°. The B y (x,t) does not show the spatial modulations with k > ki, which the 
B y (x,t) in Fig. [2] does prior to its saturation at t ps 50. Contrary to the Fig. [2H), the 
B y (x,t = 50) and the E x (x,t = 50) are obviously correlated in Fig. HH). The E x = 
whenever [B y dB y j dx) = 0, which suggests the magnetic pressure gradient as the origin 
of E x . The dominant oscillations of E x , B y in Fig. HH) are in the k 2 and k\ mode, 
respectively. The E x component evidences furthermore the presence of harmonics. 

Figure [5] illustrates the phase space distributions fi^(x,p x ) and f(x,p z ) at the 
times t = 50 and t = 120 in the simulation L 2 . The density maxima are shifted by 
L 2 /2 at both times and no further density maxima occur. The phase space structures 
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Figure 5. The 10-log of the phase space densities in units of CPs at t = 50 (a-c) 
and t — 120 (d-f) in the box L2: Panels (a,d) show the f(x,p z ) with po = m e VbT(vb). 
The beam temperature along p z is unchanged and the spatial oscillations < p z (x) >i,2 
of the beams is weak. The density oscillates by the factor w 10 2 . The fi(x,p x ) is 
shown in (b,e) and f2(x,p x ) in (c,f). Both beams separate in space and (e,f) reveal 
cool electron clouds immersed in hot electrons with a thermal spread ss Pq. 



of both beams at t = 50 reveal a high degree of symmetry, evidencing a dominant single 
filament. A second, independently developing filament would break such a symmetry as 
in the simulation L%. This is in line with the growth of B y (x, t) in Fig. H]that shows only 
a modulation with the wavenumber k\. The constant slope of ^(0.8 < x < 1.1, t = 50) 
in Fig. Hid) corresponds to a spatially uniform distribution of f(x,p x ) in Fig. [5b). The 
harmonics of E x in Fig. 0H) are related to the phase space structures at the filament 
boundaries. The f{x,p x ) at t = 120 shows a heated electron population similar to that 
in Fig. [3j The dense electron component in L2 is, however, cooler and it shows no vortex 
structures. The core populations of both beams in Fig. [5b, f) are not overlapping, as 
the current filaments do in the Fig. [3j 

The supplementary movie 2 animates in time the 10-logarithmic phase space 
distributions fi{x,p x ) and f\{x,p z ) of the beam 1 in the simulation L 2 . The formation 
of the filaments is demonstrated. The phase space evolution shows that the distribution 
f(x,p x ) contains fewer vortices and that the vortices are spread out over a smaller 
interval of p x than in the simulation L\. The spatial modulation of < p z (x) >i of the 
beam 1 is also less pronounced. The plasma thus appears to be less turbulent than 
that in L\, which may explain the more obvious relation between B y and E x in Fig. H] 
compared to the Fig. El The spatial width of the plasmon containing the dense bulk of 
the electrons of fi(x,p x ) oscillates in time. The overlap of the filaments in Fig. [5fe,f) is 
thus time dependent and related to the oscillating electrostatic field in Fig. 0b). 
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Figure 6. The field amplitudes in the box L3: The panels (a-c) show B y , E x and E z , 
respectively. The B y -amplitude reaches a steady state. The E x and E z components 
are oscillatory in space and time. The E x oscillates in space twice as fast as B y and 
both correlate well for 55 < t < 125. The E z shows a phase shift of 90° compared to 
By and its amplitude decreases in time. The B y , E x fields at t — 56 are displayed by 
(d) and we find that E x = when B y dB y /dx = 0. 



3.3. Simulation L 3 

The turbulence is reduced further, by decreasing the box length from L 2 = 2X S to L 3 = 
0.89A S . We exploit this to examine quantitatively and in more detail the relation between 
the electric and magnetic fields and their effect on the particle trajectories. Figure [6] 
displays stationary field structures. The B y (x,t) oscillates with the wavenumber k\ in 
space. The E x oscillates with the wavenumber k 2 and it is practically undamped. Both 
fields display a persistent correlation. The E z component is damped in time and it 
is shifted by the phase 90° relative to B y when the FI saturates. The damping of E z 
does not visibly influence the B y and E x , suggesting that E z is driven only during the 
growth phase of the FI and that it decouples upon saturation. The E x , B y fields show 
an excellent qualitative match between E x = and B y dB y j dx = when t = 56. 
The force of the magnetic pressure gradient on a current is expressed as 

J x B = -V B 2 /2. (4) 

The derivatives in the y and the z directions vanish in our ID geometry, B y ^> B x , B z 
and J z ^> J x , J y . The magnetic pressure gradient force on the electrons can only be 
mediated through an electric field force along x. This electric field for the normalized 
electron charge —1 is then given by Eb = —B y dB y /dx. The E x (x,t > 56) oscillates in 
Fig. [6] in space and time between E x = and an extremum. The peak field moduli 
w 0.02 are most suitable for a comparison with Eb, due to the high signal-to-noise ratio. 

Figure UK) displays the E x (x, t = 56) when the FI has just saturated and compares 
it with Eb- It turns out that E x (x) ~ 2E B (x) at t = 56, when the peak amplitude of E x 
is reached. The electric field amplitudes can be averaged over the time interval t\ = 56 
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Figure 7. (a) The E x {x,t = 56) (dashed curve) and 2Eb (solid curve), (b) The E x (x) 
(dashed line) and Eb (solid line), which have been averaged over 56 < t < 125. (c) 
The number densities for i = 56, normalized to 2n e , of both beams separately (beam 
1 is almost confined to 0.2 < x < 0.6) and both densities added together. The total 
beam density is modulated by about 30%. The density of the individual beams varies 
by an order of magnitude. 



to t 2 — 125, for which the field structures do not convect and oscillate around a constant 
mean field, to give E x = (t 2 — tx)~ / t * a E x (x, t) dt and E B = (t 2 — tx)~ / t * a E B (x, t) dt. 
The E x and Eg match in Fig. [7b). The magnetic pressure gradient is thus the origin 
of the electrostatic field. The system is oscillating around the equilibrium because 
our initial conditions did not correspond to a steady state configuration. The same 
oscillations of E x around a mean field as well as the correlation between B y and E x 
can also be observed in Fig. [2] and Fig. H] at late times for the larger boxes. The 
Eb correlates well with the E x in the simulations L\,L 2 , although the curves do not 
match as accurately as in Fig. [7J This is due to the more turbulent plasma and because 
the convection of the structures imposes either constraints on the integration time or 
requires a transformation into a moving frame, the speed of which has to be estimated. 

Figure [7b) shows the normalized number density distributions ni t2 (x) = 
(2n e ) 1 J fi,z(x, p) dp of the individual beams and also the summed distribution n\{x) + 
n 2 (x) at t — 56. The spatial oscillation period of either ra 1 (x) or n 2 (x) is L 3 , while that 
of the E x (x, t = 56) is L 3 /2. The phase space structures of electron phase space holes in 
an unmagnetized plasma would have the same periodicity as the electrostatic field [28J. 

Figure [8] shows the phase space distribution of the electrons in simulation L 3 at the 
times t = 56, 120. The < p z (x) >i^ 2 are practically unchanged along the x-direction. 
A spatial modulation is caused by the E x B-drift, the speed of which is given by the 
product of E x and B y in our ID simulations. The amplitudes of B y and of E x both 
increase according to the Figs. El [Hand [6] as the filaments get larger and, hence, the drift 
speed contribution to p z . The amplitudes of E x in the simulations Li,L 2 are several 
times the one in the simulation L3 and their drift electric fields Vf,B y are larger. The 
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Figure 8. The 10-logarithmic phase space densities in units of CPs at t = 56 (a-c) 
and t = 120 (d-f) in the box L^: Panels (a,d) show the f(x,p z ) with po = m e VbT{vb)- 
The temperature along p z and the < p z {x) >i :2 of each beam are unchanged. The 
density oscillates by the factor w 10. The fi(x,p x ) are shown in (b,e) and the f2(x,p x ) 
in (c,f). Both beams separate in space and (e,f) reveal cool electron clouds immersed 
in a tenuous electron population with a thermal spread ~ 0.4po- 

simulation boxes are also longer. The electrostatic potentials in the simulations Li,L 2 
thus exceed by far that in the simulation L 3 and the electrons can reach higher kinetic 
energies. Consequently, the spread in p x of the phase space distributions of both beams 
in the simulation L 3 is less than half of that in the simulations Li, L 2 at t — 120. 

The supplementary movie 3 animates in time the evolution of the 10-logarithmic 
distributions f\(x,p x ) and fi(x,p z ) in the simulation L 3 . The fi(x,p z ) evidences that 
the electrons are re-distributed along x, but not along p z . The electron flow along x 
oscillates. The f\(x,p x ) after the saturation of the FI has a dense electron core, which 
rotates in the o^p^-plane around x = 0.4. Two vortices in the dense electron core, 
presumably electron phase space holes, are convected with this rotating flow. 

All simulations evidence that a core of cool electrons is spatially confined. Their 
circular phase space motion in the movies around the equilibrium points x e with 
E x (x e ) = B y (x e ) = 0, for example x e = 0.4 in the movie 3, furthermore reveals, that 
they are trapped by an electrostatic potential in the x,p x plane. We consider the quasi- 
equilibrium established in the simulation L 3 after the FI has saturated. We calculate 
E x = (t 2 — ti)~ E x (x, t) dt and B y = (t 2 — t\)~ J t * 2 B y (x, t) alt for the simulation L 3 . 
We integrate both fields from t\ = 68 to t 2 = 125, when the equilibrium is established. 
The weak modulation of p z and thus v z of the CPs in L 3 allows us to estimate the drift 
electric field as E D = v^By and the total electric field E T — E x + E D . The average 
potentials Uj(x) = Uoj + Jq Ej(x)dx with the indices j = x,D,T are calculated from 
these fields and U j is set such that Uj(x e = 0.4) = 0. The potentials are expressed 
in Volts, allowing a straightforward comparison with the particle kinetic energies. The 
average fields and potentials are displayed by Fig. [91 The electric fields are such that E x 
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Figure 9. The fields Ej and the potentials Uj, both averaged over the time interval 
68 < t < 125: a) shows the electric field E x (dashed), the Ed = VbB y (dash-dotted) 
and Et = E x + Ed (solid line). Positive Ej accelerate electrons into the negative 
^-direction, b) shows the potential U x (dashed), Ud (dash-dotted) and the Ut (solid). 
The potential at x = 0.4 is the reference potential. 

destabilizes the equilibrium position x e = 0.4, because the negative E x (x > x e ) close to 
x e accelerates the electron in the positive direction and the positive E x (x < x e ) close to 
x e in the negative direction. The Ed(x ~ x e ) is confining the electrons around x ~ x e . 
The \Ejj\ > \E X \ for x m x e and Et is thus a confining force. However, the electron 
acceleration at x w x e is decreased and increased at larger \x — x e \. This is reflected also 
by the potentials. The magnetic potential invoked by Ref. [H] is dominant. However, 
the electrostatic field flattens the bottom of the potential and steepens the walls. This 
modified potential results in a bouncing time of electrons that differs from that in Ref. 
[14] . The asymmetry of the potential close to the stable equilibrium x e ~ 0.4 of the 
beam 1 and the stable equilibrium x ~ 0.85 of the beam 2 arises from the dependence 
of Ed on the beam velocity. The velocity is Vb for beam 1 and — 1>& for beam 2. The E x , 
on the other hand, acts on both beams in the same way. 

The CPs of the beam 1 should follow almost straight paths close to x e = 0.4 
and they should be rapidly reflected for \x — x e \ > 0.2. The potential difference 
Ajj = maxipT) — minlUr) ~ 1700 V in the simulation L3 should trap electrons with 
speeds up to A„ = (2eA u /m e ) 1 ' 2 /v b ~ 0.27. This matches the momentum spread of the 
cool core population in movie 3 or Fig. [SJ The oscillations of E x in Fig. [6] explain the 
periodic release of electrons from this cool core seen in the movie 3 and the oscillatory 
force imposed on the electrons by E x contributes to their heating. 

The trajectories of two CPs in the confined structure of beam 1 are followed in time 
in Fig. [10] in order to compare their dynamics to that expected from Et- The red circle 
denotes the time, when the CPs start interacting with the fields and the trajectories are 
followed until t = 125. The CP 1 has a low initial modulus of p x and CP 2 a high one. 
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Figure 10. The particle trajectories of two selected CPs: (a,c) show the x,p x and 
x,p z diagrams of the CP 1. (b,d) show the corresponding diagrams for the CP 2. 
The red circle denotes the starting point of the trajectory. Both CPs follow straight 
paths in the x,p x plane for 0.33 < x < 0.47 and they are rapidly reflected outside this 
interval. The E x B y drift imposes diagonal paths of the CP 2 in the x,p z plane. 



Both CPs follow straight paths in the interval 0.33 < x < 0.47, in which Ej- in Fig. [9] 
is small. The phase space path of the faster CP 2 is smoother than that of CP 1. The 
low speed of CP 1 implies a long crossing time of the interval with a low modulus of E? 
and the CP 1 experiences several oscillation cycles of E x . The simultaneous action after 
the saturation of the FI of the quasi-stationary B y and the oscillatory E x , which both 
vary in space, implies that the electron acceleration is a function of the position and of 
the phase of E x relative to B y . This phase- dependence results in a randomization of the 
particle trajectories, contributing to the plasma heating. The faster CP 2 crosses the 
bottom of the potential more rapidly and it experiences lower relative speed changes by 
the E T . The particles are reflected outside the interval 0.33 < x < 0.47 by the \E T \. 
Both CPs are trapped because their speed is less than that required to overcome A[/. 

The i? 2 -fields are weak but not negligible. A comparison of the Figs. |2j HI and [6] 
demonstrates, that its time-evolution depends on the box length. The E z grows in all 
simulations prior to the saturation of the FI and its phase is shifted by 90° relative to B y . 
It is thus pumped through Ampere's law by the growth of J z . The initial, low-frequency 
oscillations of E z with k = k\ damp out in all simulations. They are replaced in the 
simulations L\, L 2 by faster oscillations that have k = and k = k\. This is not the case 
for the E z in the simulation L3. The power spectrum of E z reveals the origins of the 
damping of the low-frequency modes and of the growth of the high-frequency modes. 
This power spectrum is obtained by Fourier transforming E z (x,t) — > E z (k,uj) over the 
full box and simulation time. The power spectrum E(k,u) = \E z (k,uj)\ 2 , which we plot 
for the simulations L 2 , L 3 in Fig. [TUJ 

The E(k, oS) reveals the following for the simulation L3. The spectrum of E z consists 
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Figure 11. The 10-logarithmic power spectra of E z , integrated over < t < 125 and 
normalized to their peak value. Frequencies are normalized to fl p . (a,c) correspond to 
the simulation L% and the k is normalized to 27r / L3. (b,d) correspond to the simulation 
L2 and k is normalized to 2n /L2. (a) reveals a discrete wave spectrum in k while the 
spectrum in (b) is turbulent. The frequency spread is due to the growing amplitude 
of E z . The frequency of the fast X-mode with lu ck\ is indicated by vertical lines in 
(c,d). Only (d) shows a pumping of the fast X-mode. 



of the fundamental mode at k — k\ and of secondary waves at kj with j=l+2n and 
n = 1,2, . . The modes with the k > k\ are not harmonics of k\. They can thus not 
originate from a self- interaction of E z or from a coupling of E z with B y , which has 
also k — k\. It appears that E z is interacting with the k 2 mode of E x and, possibly, 
also with By. The power spectrum is spread out in frequency, because E z initially 
grows exponentially and is thus not monochromatic. The E z can couple to the fast 
extraordinary (X-)mode with the same polarization. Its frequency u ~ ck in physical 
units at k = ki is high for the simulation L3 due to the high k\ of this simulation and 
no power is transfered from the low-frequency turbulence to the high-frequency fast 
X-mode. The transient oscillations of E z in the simulation L3 damp out. 

As we go from L3 to L2 the frequency of the fast X-mode at k — k\ is lowered by a 
factor 2.25. Figure [11] evidences that now wave energy is coupled into the fast X-mode 
at k = k\. This pumping also transfers wave power to the mode with k = 0. This may 
be achieved by a parametric interaction of the low-frequency mode at k = ki and the 
high-frequency mode at k — k\. The frequency spread of the initial modes implies, that 
it is not necessary to fulfill exactly the frequency and wavenumber conservation of the 
3-wave interaction. The fast X-modes are linearly undamped, explaining the persistent 
high-frequency oscillations of E z in the Figs. [2] and HI 
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4. Discussion 

In this work we have examined in detail the saturation of the electron beam filamentation 
instability (FI) with one- dimensional PIC simulations. The one-dimensional geometry 
implies that we can not model the FI beyond its saturation, because the filament merging 
is excluded [16] . However, we can readily distinguish between the electrostatic fields with 
their polarization vector along the simulation direction and the electromagnetic waves. 
The phase space distribution can also be sampled with a high statistical accuracy. We 
have modelled here only mobile electrons. The simulation compensates their charge by 
introducing an immobile positive charge background that cancels the electron charge. 

The FI grows magnetic fields out of noise by a spatial separation of the currents 
of the, initially uniform, electron beams. A broad wavenumber spectrum is destabilized 
and the amplitudes of the individual filamentation modes grow aperiodically and 
exponentially. The magnetic field amplitude can not grow indefinitely. The magnetic 
amplitude saturates, once the bouncing frequency of the electrons in the magnetic 
potential becomes comparable to the growth rate of the instability [H] . More recently, it 
has been pointed out that the electrostatic fields upon saturation may not be negligible 
[IB] . These electrostatic fields have been connected qualitatively to the magnetic 
pressure gradient [20], [21] . The observation, that the saturation of the FI is qualitatively 
unchanged by a beam-aligned magnetic field, has led to the suggestion in Ref. [20], that 
it is not the magnetic trapping but the electrostatic fields that quench the instability. 
This is, because the bouncing frequency of electrons depends on the strength of the 
flow-aligned magnetic field, while the electrostatic field does not; the pressure gradient 
of a uniform magnetic field vanishes and leaves unchanged the electrostatic field. 

The electron beams in our simulations have been equally dense and cool. The 
beams counter-propagated at a non-relativistic speed orthogonally to the simulation 
direction. Picking a box length comparable to the filament size has allowed us to obtain 
a quasi-monochromatic wave spectrum and one pair of filaments. We have demonstrated 
quantitatively, that the electrostatic field is indeed driven by the magnetic pressure 
gradient. The initial conditions we have used here do not represent a steady state 
and the saturated plasma oscillates around its equilibrium state. The amplitude of the 
electrostatic field oscillates around the mean value, which is expected from the magnetic 
pressure gradient. We have demonstrated this for the small filament. Small means here, 
that no sub-structures can form because thermal effects [HI [29] limit the maximum 
unstable wavenumber. We have quantified the relative importance of the electrostatic 
field and of the magnetic field. The electrostatic field is comparable in strength to the 
drift electric field v\,B y and both are thus relevant for the saturation of the FI and for the 
selected plasma parameters. The electrostatic field oscillates in space twice as fast as the 
magnetic one. Its effect is to reduce the force on the electrons close to the equilibrium 
point of the respective filament and to increase it farther away from it. The effective 
potential obtained from a summation of the magnetic pressure-driven electrostatic field 
and of the drift electric field is thus not a cosine. 
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The electrons can move almost freely through the potential well and are reflected 
at a relatively thin layer, as we have demonstrated for two representative computational 
particles. The effective potential is not symmetric with respect to both filaments. 
This is, because the drift electric field depends on the beam flow direction, while the 
magnetic pressure-driven electric field acts on all electrons the same way. The wall 
of the potential that is confining one filament is located well inside the potential well 
of the second filament. Filaments thus remain only separated if they have opposite 
beam flow directions, through which they can overlap without merging. This effective 
potential will also accelerate ions [18] . This ion acceleration is, however, exaggerated 
in ID simulations where the potential is stationary in space. The filaments merge in 
higher-dimensional simulations and the potentials are not longer stationary. We have 
thus deliberately omitted to consider ion effects. 

The electrons are heated up when the potential grows and also by their simultaneous 
interaction with the oscillatory electric field and the steady magnetic field. A dense and 
spatially confined electron population (plasmon) maintains the current responsible for 
the magnetic field after the heating has taken place. A hot electron population overcomes 
the potential well and presumably contributes to the quenching of the FI, which can be 
suppressed by a high temperature orthogonal to the beam flow direction [29] . 

We have then assessed the impact of the filament size on its dynamics. The 
probability distribution for the filament size has been sampled using a long one- 
dimensional simulation box in Ref. [21]. We have examined here two other filaments, 
the size of which we have selected according to the probability distribution. The size 
of the medium filament is close to the most common size. The large filament has a size 
close to the maximum observed one. Larger filaments can only form, if filaments merge 
in higher-dimensional systems [TH [23] because the growth rate of the FI decreases, as 
we go to smaller wavenumbers. Long waves are then outgrown by shorter waves. The 
electrostatic fields of the medium and the large filament have been correlated with the 
magnetic field after the FI has saturated and they have oscillated around an equilibrium 
value. This equilibrium value is close to the magnetic pressure gradient- driven field, as 
for the small filament. We have not shown it, because the fit is not as accurate as for the 
small filament, due to the drift of the structure and the increased levels of turbulence. 
The sub- structure, i.e. the merging of smaller filaments in the large filament, prior to 
the saturation of the FI has furthermore prevented a clear correlation of the electrostatic 
with the magnetic field at the saturation time. Both the medium and large filaments 
propagated after they have saturated, even though the phase speed of the waves driven 
by the FI is zero. The speed has, however, been less than the initial beam thermal speed. 
A filament can be accelerated by the saturation of the FI. Any net momentum of the 
hot and untrapped electron population must be balanced by an oppositely directed net 
momentum of the trapped electrons. The electromagnetic fields are tied to the trapped 
electrons and convect with the plasmons. The direction of the convection is random. 
The movie in Ref. |21j shows that different convection speeds of neighboring filaments 
result in their spatial oscillation rather than in a convection. 



Saturation of the filamentation instability 



18 



The mean momentum along the beam flow direction has been spatially modulated. 
Our comparative study of three filament sizes has shown that the magnetic field 
amplitude and the electrostatic field amplitude both increase with the filament size. 
They modulate the mean velocity of the beam through the E x B-force. This modulation 
thus increases with the filament size, explaining the observation in Ref. [23] that trapped 
electrons reach increasingly higher speeds as the filaments merge. 

The magnetic field has varied linearly in space over wide spatial intervals in all case 
studies. A magnetic field amplitude with a constant gradient results also in a magnetic 
pressure-driven electric field with a constant gradient. The Fourier transform over space 
of a curve with a constant gradient results in a power-law spectrum. This is observed 
also in longer- and in higher- dimensional simulation boxes [211 1231 130] . 

We have examined the electric field component along the beam flow direction. Its 
wavenumber, the phase shift relative to the magnetic field and that it grows during 
the linear phase of the FI implies that it is driven through Ampere's law. These 
low-frequency oscillations are transient modes and they damp out. However, we have 
confirmed here a previous observation [18], that the fast extraordinary mode is pumped 
by this wave component. We could observe this only for the two large filaments. The 
finite box size introduces a discrete wave spectrum. The larger the box length, the lower 
the frequency of the electromagnetic mode. The frequency has been low enough in the 
large simulation boxes to absorb wave energy from the low-frequency turbulence. The 
discrete wave spectrum is, however, a finite box (numerical) effect. 
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